function [schidx,nsidx] = TwoChoice_StimDelayTraining(sch,schidx)
% [schidx,nsidx] = TwoChoice_StimDelayTraining(sch,schidx)
%
%  Extend stimulus delay (cue_delay) and response window by adaptive
%  procedure.
%  
%  If correct response, increment cue_delay by DELAY_STEP, otherwise
%  decrement cue_delay by DELAY_STEP.  Delay boundaries are set by
%  DELAY_MIN and DELAY_MAX. Initial delay is set by DELAY_START and the
%  number of correct responses in a row to increment cue_delay is set by
%  DELAY_ADV.
%
%  Use the following parameters in schedule.writeparams including *
%  ('Write' or 'Read/Write' functionality in ScheduleMaker):
% 

%  Takes the schedule and schidx as input.
%       SCH is the schedule structure
%       SCHIDX is an Nx1 matrix with the number of rows equal to the number
%           of unique trials; i.e. size(trials,1).
%  Returns an updated version of SCHIDX and NSIDX
%       NSIDX is the next schedule trial index, i.e. the next index
%           selected from the SCH.trials matrix to run.
%       SCHIDX should be updated with the number of runs for each unique
%           trial in SCH.trials.  The returned version should have
%           SCHIDX(NSIDX) incremented.
%
%  This same format for input and output paramters must be met for custom
%  trial select functions (i.e., experiment.trialselectfcn);
%
%  DJS (c) 2011

global G_RP 
persistent GSEQ GSEQseed TCNT


if isempty(GSEQ)
    % Generate Gellermann (1933) randomized sequence of trials
    %% Gellermann, J1933
    n = 10;
    p = zeros(2^14,n);
    for i = 1:n
        p(:,i) = randperm(2^14);
    end
    
    m = [-ones(2^12,n); ones(2^12,n)];
    m = unique(m(p),'rows'); clear p
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    %% rule 1
    ind = sum(m,2) ~= 0;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    %% rule 2
    tm1 = m > 0;
    tm2 = m < 0;
    for i = 1:n-4
        ind = sum(tm1(:,i:i+3),2) > 3;
        ind = ind | sum(tm2(:,i:i+3),2) > 3;
        tm1(ind,:) = [];
        tm2(ind,:) = [];
        m(ind,:)   = [];
    end
    
    %
    ind = false(size(m,1),n-2);
    tm1 = m > 0;
    tm2 = m < 0;
    for i = 1:n-2
        ind(:,i) = sum(tm1(:,i:i+2),2) == 3;
        ind(:,i) = ind(:,i) | sum(tm2(:,i:i+2),2) == 3;
    end
    ind = sum(ind,2) >= 2;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    %% rule 3
    ind = abs(sum(m(:,1:n/2),2))>1;
    ind = ind | abs(sum(m(:,n/2+1:end),2))>1;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    %% rule 4
    dm = diff(m,1,2)/2;
    ind = sum(abs(dm),2) ~= n/2;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    %% rule 5
    % rule 5 is implied
    
    %% Construct trial sequence
    a = m(1:size(m,1)/2,:);
    b = m(size(m,1)/2+1:end,:);
    
    ridx = randperm(size(a,1)); a = a(ridx,:);
    ridx = randperm(size(b,1)); b = b(ridx,:);
    
    m = [a fliplr(b)];
    
    GSEQ = reshape(m',1,numel(m));

    % Start from random position within sequence
    % Since Gellermann rules for randomization only allow for 44 sequences
    % of 10 trials yield 440 trials total, limit seed to first 100 trials
    % in the sequence.
    GSEQseed = randperm(100);
    GSEQseed = GSEQseed(1);
end


boxid = sch.boxid;

% find some paramters
% sdidx = findincell(strfind(sch.readparams,'cue_delay'))+1;
rvidx = findincell(strfind(sch.readparams,sprintf('response_code~%d',boxid)))+1;

% stimulus delay parameters
dstartidx = findincell(strfind(sch.writeparams,sprintf('DELAY_START~%d',boxid)));
dstepidx  = findincell(strfind(sch.writeparams,sprintf('DELAY_STEP~%d',boxid)));
dadvidx   = findincell(strfind(sch.writeparams,sprintf('DELAY_ADV~%d',boxid)));
dminidx   = findincell(strfind(sch.writeparams,sprintf('DELAY_MIN~%d',boxid)));
dmaxidx   = findincell(strfind(sch.writeparams,sprintf('DELAY_MAX~%d',boxid)));

DELAY_START = sch.trials{1,dstartidx};
DELAY_STEP  = sch.trials{1,dstepidx};
DELAY_ADV   = sch.trials{1,dadvidx};
DELAY_MIN   = sch.trials{1,dminidx};
DELAY_MAX   = sch.trials{1,dmaxidx};

RV = sch.response_vals;

trials = sch.trials;

% which module?
writemod = sch.writemodule(1);
RP = G_RP(writemod);

% correct side: -1 Left; 1 Right
csidx = findincell(strfind(sch.writeparams,'CorrSide')); 
corrside = cell2mat(sch.trials(:,csidx));

if isempty(RV)
    tidx = 0;
    TCNT(boxid) = 1;
else
    tidx = size(RV,1);
end

subcorrside = find(corrside == GSEQ(GSEQseed + tidx));

% initialize first trial
if isempty(schidx)
    schidx = zeros(size(trials,1),1); % initialize SCHIDX
    SD = DELAY_START;
else
    SD = RP.GetTagVal(sprintf('cue_delay~%d', boxid)); % last stimulus delay
end

% last N responses during response window (light on)
if size(RV,1) >= DELAY_ADV && TCNT(boxid) >= DELAY_ADV
    HITind = bitget(RV(end-DELAY_ADV+1:end,rvidx),6);
    HR = sum(HITind)/DELAY_ADV;
    if HR >= 0.8
        SD = SD + DELAY_STEP;
        
    elseif HR <= 0.2
        SD = SD - DELAY_STEP;
    end
    
    TCNT(boxid) = 1;
else
    TCNT(boxid) = TCNT(boxid) + 1;
end

% respect bounds
if SD > DELAY_MAX, SD = DELAY_MAX; end
if SD < DELAY_MIN, SD = DELAY_MIN; end

% Update RPvds Tags
RP.SetTagVal(sprintf('cue_delay~%d', boxid),SD);

if TCNT(boxid) == 1
    fprintf('BOX %d: STIM DELAY = %1.0f\n',boxid,SD)
end

% update SCHIDX
% give priority to least chosen trials
i = min(schidx(subcorrside));
i = find(schidx(subcorrside) == i);
r = randperm(length(i));
nsidx = subcorrside(i(r(1)));
schidx(nsidx) = schidx(nsidx) + 1;









